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We show how the worldline quantum Monte Carlo procedure, which usually relies on an artificial 
time discretization, can be formulated directly in continuous time, rendering the scheme exact. 
For an arbitrary system with discrete Hilbert space, none of the configuration update procedures 
contain small parameters. We find that the most effective update strategy involves the motion of 
worldline discontinuities (both in space and time), i.e., the evaluation of the Green's function. Being 
based on local updates only, our method nevertheless allows one to work with the grand canonical 
ensemble and non-zero winding numbers, and to calculate any dynamic correlation function as easily 
as expectation values of, e.g., total energy. The principles found for the update in continuous time 
generalize to any continuous variables in the space of discrete virtual transitions, and in principle 
also make it possible to simulate continuous systems exactly. 
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I. INTRODUCTION 



(N : 

■ Quantum Monte Carlo (MC) simulation is the most powerful available method, if not the only one, of obtaining 
accurate results for complex systems, where analytic solutions are not possible and exact diagonalization methods 
do not work because of the enormous Hilbert space. However, most MC schemes are. far from ideal, and suffer from 
' significant shortcomings. These include (see, e.g., the most recent review article Ref.El) 

a) systematic errors due to artificial time discretization, which in most schemes scales as (At) 2 , where At is the 
^ — , ' time slice width; 

Q\ , b) restriction of the simulation to the zero winding number sector M — (a configuration in which world lines 
connect the initial state |ai,aa, ■ ■ ■ ,ai) at t = to the final state |7i,72, ■ ■ ■ ,7l) at t = /3, with the set {7^} being 
obtained by cyclically permuting {a^} M times (and all topologically equivalent configurations), is said to have a 
winding number M). Such a restriction results in systematic errors too, which however vanish with increasing system 
size. Also, one loses the ability to study topological excitations in the system, e.g., vortices or supercurrent states; 

c) working with a fixed number of particles N — const (canonical ensemble); 

d) critical slowing down problem, which arises close to a second-order phase transition. This problem is closely 
, related to constraints (b) and (c), and is indicative of inefficient procedures used to update configurations with large 

length scales; 

e) slow accumulation of statistics when calculating correlation functions of operators not present in the initial 
Hamiltonian, e.g., the Green's function; 

f) small acceptance rates in update procedures. These may be due to from small parameters present in the 
formulation of the MC scheme, or systems described by Hamiltonians with different energy scales (e.g., when the 
hopping matrix element t is much smaller than the typical potential energy change U 3> i) , or the necessity of global 
Metropolis updates, which arise in certain cluster-update algorithms. 

g) anomalous dependence of the computation time on system size (due to self-averaging effects in the thermodynamic 
limit, the computation time required to achieve given accuracy is expected to be system-size independent); 

h) notorious sign problem, which emerges when the configuration weight is not positive definite. Since we do not 
see any reasonable solution of the sign problem in the general case, in what follows we exclude it from the discussion. 

To eliminate some of these shortcomings, a number of different MC schemes were developed. Unfortunately, none 
of the existing schemes succeeded in solving all of them (leaving the sign problem aside) in the general case: there 
are extremely efficient algorithms which are far from universal, while the efficiency of existing universal algorithms is 
far from high for a large number of problems. 

The standard worldline algorithm is based on imaginary time discretization and utilizes the small parameter t At <C 
1 in an approximate treatment of noncommuting operators in the Hamiltonian, known as Trotter break-upBu. Physical 
intuitiveness and easy programming probably make this method the one most widely used. On another hand, its weak 
points range over the whole list from (a) to (f), the most severe ones being (e) and (f). 

In the worldline algorithm, one describes the configuration by specifying the system state \ak) at all time slices 
Tk = k At, where ft = 0, 1, ... , Kp and txg — = (3- The system state is then conventionally defined in the basis 
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set in which the potential energy of the system is diagonal, i.e., in the site representation. Consider, as a typical 
example, the Hamiltonian of interacting particles on a lattice 

H = -t 4 a j + UijUiUj , (1.1) 

<ij> ij 

where a\ creates a particle on site i, t is the hopping matrix element, ni = ct|<Xj, and < ij > denotes nearest- neighbor 
sites. From now on, we call points in time at which the system changes state "kinks." The typical separation in 
time between two adjacent kinks on the same site is of order 1/t and independent of At, so that for small At there 
are some l/(iA-r) S> 1 time intervals between them. The acceptance rate of the variation suggesting creation of a 
new kink - antikink pair is proportional to the square of a small parameter, (£At) 2 . On the other hand, when the 
MC procedure suggests shifting an already existing kink to the nearest point in time, the corresponding variation of 
the configuration is accepted with probability ~ 0(1). Thus, on average, by selecting different random time slices, 
it takes some 1/(£At) 2 attempts to create a new kink - antikink pair and 1/(£At) attempts to move a kink to 
the nearest position in time. Still, the updated configuration is only slightly different from the previous one, and 
expectation values calculated before and after the variation are strongly correlated. An uncorrelated contribution 
of the given configuration fragment is obtained by shifting the kink a distance of order 1/t, which requires some 
l/(tAr) 3 operations, since the kink shift process is diffusive in nature. This means that the autocorrelation time in 
the standard worldline algorithm grows as oc (At) -3 even in the absence of critical slowing down. Since all update 
procedures are local, the algorithm is subject to critical slowing down near the transition temperature. 

In order to calculate the Green's function Q(i, r), two worldline discontinuities are inserted at time slices t\ = and 
T2 = t (in other words, one extra worldline is inserted in or removed from the interval [ti,T2] )au. One then probes 
different configurations using standard update procedures and collects statistics in a d-dimensional histogram, which 
describes the spatial separation i between the discontinuities. The length of the time interval is then changed, and 
the same calculation is repeated. One MC step, i.e., the number of update operations performed between successive 
"measurements," whereupon another point is included in the statistics of the calculated quantity, is proportional to 
L d f3, where L d is the number of lattice sites considered. Thus it requires about L d (3 operations to include only one 
point in the (d + l)-dimensional spacetime histogram for G. 

In a way, the standard worldline procedure of calculating Q(i,r) has an anomalous dependence on N and /?, since 
it takes at least (N(3) 2 operations to update the whole histogram (typically G decays in space and time, and large 
scale behavior requires much more computation). We note, that winding numbers and the grand canonical ensemble 
average can be incorporated, in principle, in the worldline algorithm. It is sufficient to consider separate contributions 
to the statistics of Q(i, r) when i = ML and r = n(3 with integer M and n. However, in practice, only small systems 
at rather high temperatures can be considered using this algorithm. _ _ 

The determinant method based on the Hubbard - Stratonovich transformationDcl also uses the discrete-time Trotter 
break-up, and thus becomes more and more inefficient due to long autocorrelation times when At — ► (points (a) 
and (f) above). It has an important advantage over the worldline method in calculating the Green's function, since 
it works with the grand canonical ensemble. However with increasing system size, the calculation time scales as L 3 
(point (g)) and some of the procedures become ill-conditioned at low temperatures. 

Another technique allowing Green's function calculations is called Green-function MC (or, more generally, the 
projection-operator method)0. It is applicable at zero temperature only, and the final result for Q{i,T = 0) depends 
on the trial wave function (we are not aware of whether it is possible to calculate the time dependence of Q{i,t) by 
this method). nnn n 

The stochastic series expansion (SSE) techniqueEra'E2l, which stems from the Handscomb's methocO, relies on the 
direct Taylor expansion of the statistical operator. This scheme is exact (contains no systematic errors). SSE has 
clearly demonstrated that time discretization is an artificial trick that is not at all necessary for MC simulation. Since 
an elementary update in the SSE scheme is equivalent to roughly l/(tAr) 3 updates in the standard worldline method, 
it results in a significant drop in computation time for high-precision calculations. The rest of the problems, i.e., (b)-(f), 
survive in the SSE approach (point (f ) still applies, because by expanding in powers of the full Hamiltonian one h as t o 
compare weights corresponding to the kinetic- and potential-energy terms, and if, e.g., U 3> t in the Hamiltonian (tL.lt) , 
then small acceptance rates appear in the update procedures). Still, away from the transition point, for large systems 
at low temperature, for which U ~ t, the SSE method is superior in evaluating basic thermodynamic properties like 
the total energy and density - density or current - current correlation functions. ._. ._. 

A qualitatively new class of extremely efficient MC schemealj has been developed in recent yearsEftllj. These 
schemes are based on the so-called loop cluster update (LCU) algorithm, which performs nonlocal updates for worldline 
loops with sizes as large as the system itself. Apart from solving the problem of critical slowing down, it also allows one 
to work in the grand canonical ensemble and with nonzero winding numbers. From this method we learn that problems 
(b), (c) and (d) can be circumvented. Unfortunately, the LCU algorithm, as far as we know, is not universal. It applies 
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to spin systems and to hard-core Hubbard models, but was never formulated for the general lattice Hamiltonian, like 
interacting soft-core bosons with arbitrary f/y, arbitrary density (chemical potential), or on-site disorder. Another 
shortcoming of the LCU, which in a sense can also be called nonuniversality, is that it does not admit of a universal 
codeEHl. It should be noted that LCU allows for considerable generalization to the cases of external magnetic field 
and disorder but generally speaking the cost is lost efficiency because of the exponentially small acceptance rates for 
large loopsBl. |—| 

It is shown in Ref.113 how to build a path integral in continuous time for quantum systems in a discrete basis. The 
configuration is specified by transition times and system states before and after the transition. Within this description 
one can formally think about taking the limit At — > in the standard approach. However, to implement this 
description one has to formulate the update process. In the standard worldline algorithm, one of the basic procedures 
is generation of new kink-antikink pairs when system evolution on a given site changes from |a) — > \a) — > |a) to 
\ a ) — ^ |"y 7^ ck) — >■ \&). In the continuous limit, the acceptance rate for such a variation vanishes as (At) 2 , and thus 
the problem of qualitatively new update principles arises. |— ■ nn 

Recently, two independent continuous time schemes utilizing the ideas of RefJl3 were developedliSEH Beard and 
WieseE^I find that with the LCU algorithm, one can go directly to the continuous-time limit Ar — * 0, thus rendering 
the LCU algorithm exact. _ 

The general solution to the problem of configuration update in continuous time is found in Ref.E 2 ]. The resulting 
continuous-time worldline (CTWL) method is exact, and like SSE, is I/(tAr) 3 times more efficient than finite-Ar 
local schemes. It also completely eliminates problem (f), since none of the procedures relies on small parameters 
(potential energy is accounted for in the exponent, and one does not have to weight relative contributions to the 
statistics of the potential and kinetic energy terms, as happens in SSE). 

In its original formulation, the CTWL approach did not solve problems (b) - (e), and was tested only on a simple 
single-particle Hamiltoniar£3. In the present article, we present a complete description of the CTWL approach to the 
statistics of arbitrary many-particle system with discrete Hilbert space. We demonstrate that it enables one to solve 
problem (e) in a physically intuitive way by formulating the local update procedures in terms of the motion of two 
worldline discontinuities (in what follows we call them "worms") in space and time, i.e., in terms of a calculation of 
the Green's function. During one MC step (consisting of N/3/t operations) the whole histogram for Q(i, r) is updated, 
which means that G is calculated as efficiently as, say, the total energy, and is not affected by point (g). Since 
Q(i = [M X L X , M y Ly, . . .], n(3) with integer M x , M y , . . . and n describes a system with n extra particles and winding 
numbers {M}, we are working in the grand canonical ensemble. This solves problems (b) and (c). 

Closer examination of the loop building rulealj'E^I for the Heisenberg Hamiltonian shows their remarkable similarity 
to the evolution of an extra worldline segment. The crucial difference is that only closed loops are considered by the 
LCU algorithm, while our scheme considers all the intermediate configurations as well, and utilizes them for the Green's 
function calculation. .Working in the extended configuration space, which includes discontinuous worldlines, we use 
local Metropolis-typea updates only. However, when discontinuities annihilate, and we return to the configuration 
space of closed worldlines, the net result of the update is of global character. Since CTWL with "worm" updates 
effectively mimics single-loop LCU, we may hope that it possesses all the remarkable features inherent in LCU, and 
in particular that it solves, or at least softens, the problem (d). 

To summarize, we propose a method which is exact, complete (allows calculation of any correlation function), and 
universal (applies to arbitrary quantum systems with discrete Hilbert space, and enables one to write a unified code, 
that is simultaneously applicable to lattice bosons and arbitrary spins, with arbitrarily long-range interactions and 
disorder). The sign problem now becomes the only stumbling block on the way of making quantum MC an "ideal" 
computational tool for studying complex systems. 

This paper is organized as follows. In Section || we formulate the general principles of the continuous-time worldline 
approach. In Section IH we introduce the update procedures that we find to be the most effective and sufficient for 
simulation of the quantum statistics of many-particle systems. In Section |y| we demonstrate the advantages of the 
new method by presenting some of results that cannot be obtained by any other MC approach: the Green's function 
and the critical index of the ID boson Hubbard model at the quantum critical point, and the low-energy properties 
of the strongly-disordered Bose glass phase. In Section ^ (and Appendix A) we discuss the feasibility of increasing 
the efficiency of our method in the case of long-range interactions, and consider the feasibility of generalizing of the 
method to continuous systems. 



II. GENERAL PRINCIPLES 



Let Hq and V be the diagonal and off-diagonal parts of the Hamiltonian H in a chosen representation corresponding 
to the full set {a} of eigenstates of Hq, with Hq | a) — E a \ a). The statistical operator can then ordinarily be 
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related to the Matsubara evolution operator a in the interaction picture, i.e., we write e @ H = e P H °a, with 

<7=1- f dTV{r) + ... + {-l) m [ dr m --- [ 2 dT 1 V(T m )---V(r 1 ) + ... , (2.1) 
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where V(t) = e TH "Ve~ rHa . Without loss of generality and in accordance with typical forms of Hamiltonians of 
interest, V can be written as a sum of elementary terms Q s , whose action on any function from the set {a} results in 
another function from this set: 

V = ^2Q S , Q.|a)=-g 7a (a)| 7 ) ( 7 = 7 ( s ,«)). (2.2) 



Since V is Hermitian, for any s in the sum (2^2) there exists an s' such that Q s / = Q\. We rewrite Eq. (2.1) in 
components (below E ai = E a — Ey): 
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dr m --- / dn q a v{s m )e TmE - ■ ■ ■ qxj(si)e TlE ^ +... . (2.3) 



£ 

L,...,S» 

Note that there is no additional summation over the indices of the intermediate complete sets (labeled by Greek 
letters), since these are defined in a unique way by configurations of (si, S2, ■ ■ ■ , s m ). 

We confine ourselves to the case of finite-range interaction, which is defined by the requirement that for each term 
si of elementary operators {Q s } there exists only a finite number of terms S2 for which the condition 

[QsM, Q S2 (r 2 )] =0 (2.4) 



is not met. In the case of finite-range interaction, the structure of the series (2.3) is dr astic ally simplified, the 



simplification being of crucial importance for practical realization of our algorithm. From ( |2.4| ) it follows that up 



to an irrelevant change in the indexing of energies and matrix elements, one can ignore the chronologic al o rder of 



Q Si (ti) and Q S2 (t2) in the evolution operator. This suggests representing a general term of the series (2.3) in the 
following form. First, we introduce the notion of a "kink of type s," which is characterized by a time r, a matrix 
element q Q7 (s), and a diagonal-energy difference E ai . The former two we will refer to as parameters of the kink. It 
is essential that (i) to obtain parameters of a kink one need not know explicit ly t he whole state | a ) , or | 7 ) — local 



information is enough; (ii) to specify a particular structure of a term in Eq. (2.3), including the chronological order 
of all noncommuting operators, it is enough to specify for each kink its associated neighbors, i.e., the noncommuting 
kinks nearest in time. 



Now our goal is to describe in general terms a stochastic process that directly evaluates Eq. (2J5). For simplicity, we 
assume that all q a p (s) are positive real numbers. (In many particular alternative cases, a straightforw ard generalization 
is possible, but usually at the expense of convergence.) Summations and integrations in Eq. (|2.3| ) then can be 
regarded, up to a normalizing factor, as an averaging over the statistics of different configurations of kinks, each 
configuration being defined by a certain number of kinks of certain types, their associations and particular positions 
in imaginary time. The Monte Carlo process should examine these statistics by generating different kink configurations 
in accordance with their weights. The global process will consist of a number of elementary subprocesses, each being 
responsible for certain modifications of a particular type. 

An update procedure of a general type should involve subprocesses of creation and annihilation of kinks. Clearly, 
the qualitative difference between discrete- and continuous-time QMC schemes is associated with processes of just 
this kind. To introduce the general principles of construction of subprocesses that change the total number of kinks, 
we consider some particular (but still rather rich) class of elementary transformations (which seems to be sufficient 
for all practical purposes). By an elementary transformation we mean a subprocess which either only creates or only 
annihilates a certain number of kinks. The set of elementary subprocesses can be decomposed into self-balanced 
creation - annihilation pairs. Our task then is to specify the structure of creation and annihilation subprocesses, and 
to derive the balance equation that would guarantee that the statistics generated by each pair of subprocesses does 
really correspond to that introduced by Eq. (^) . 

Let some subprocess create n kinks of given type s%, s%, . . . , s n , the temporal positions of the kinks being specified 
by the n-dimensional vector f = {n, T2, . . . ,T n }. In the most general case, the creation procedure involves two steps. 

First, one suggests to create n new kinks at f £ T, where T is a certain region in the n-dimensional space of times 
Ti, T2, ■ ■ ■ , r n . The probability density W(f) of choosing a given r is, generally speaking, arbitrary, provided W(f) is 
nonzero at every physically meaningful configuration of kinks. 
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In the second step, one either accepts (with probability -P acc (T)) or rejects the suggested modification. 
The annihilation procedure is much simpler. The n kinks of given type Si, S2, ■ ■ ■ , s n and with r € T are either 
removed (with probability Pi em (f)) or remain untouched. 

The equation of balance for the given pair of subprocesses reads 

A Pc W(f) P acc {f)df - dA„(f)p a P rom (r) = . (2.5) 

Here Aq (A n (f)) is the probability (probability density) of finding a configuration without the specified n kinks (with 
the specified n kinks at the given f) . We have also introduced the probabilities p c and p a of addressing the creation and 
annihilation subprocedures. In the next section we will see how it can turn out quite naturally that these probabilities 
do not coincide. 



The statistical interpretation of Eq. (2.3) implies 



dMf) = drf[ q ( Sj )e^, (2.6) 



3=1 



where q(sj) = qatpAsj) and AEj = E a . — Ep-. Combining ( |2 . 5[ ) and (2.6) we obtain the necessary and sufficient 
condition for the pair of subprocesses to be self-balanced: 

W(f)P acc (f) = ^ R ^ = Y[q( Sj )e AE ^ . (2.7) 



^rcm \ 7j Pc 



Given W(f), the condition (2.7) is satisfied, e.g., by the following obvious choice of P acc and P r( 



_ . R{f)/W{f), UR{t)<W(t) , 

^ cc[T> ~ I 1 , otherwise , 1 8j 



P _ j W(t)/R(t), if R(t)>W(t) , 

^remirj - < otherwise . 1 J > 



From (2.8) it can be seen that there is a certain reason for choosing W(t) oc R(t), as in this case P acc becomes 
independent of r, and the accept - reject decision can be made before suggesting a particular configuration, thus 
saving computational time. However, if the structure of the function R(f) is complicated, the numerical generation 
of the corresponding distribution will be very expensive. In this case it is better to take W(f) oc R(f), where R(t) is 
some "coarse-grained" approximation to R{f) with a simple form. 

We do not consider here a general theory of subprocesses that do not change the number of kinks, since it is basically 
the well-known theory of taking multidimensional integrals by standard Monte Carlo procedures. Particular examples 



of such subprocesses can be found in Ref.Ej and in the next section. 



The foregoing approach does not involve any explicit truncation of the series (2.3). One might wonder, however, 
what the effect of implicit truncation in the practical realization of the process would be, due to the finite size of the 
computer memory. To this end we note that even for simulations of many-particle systems, where the typical number 



of kinks Nkink (that is, the typical number of terms in the series (2.3)) that contribute to the final result is really large, 
and one might expect the memory/ accuracy problem, the effect can be easily made absolutely negligible. Indeed, 
from the Central Limit Theorem, it follows that the number of kinks in significant configurations has a Gaussian 
distribution with the peak at Nkink and a half- width of order y/ Nkink (cf. RefM). If one just reserves at least twice as 
much memory as necessary to describe the configuration with N^n^ elements, then during a computation spanning the 
age of the Universe, the system will not fluctuate to states which cannot be fit into memory. The implicit truncation 
error thus can be made astronomically small. 

III. UPDATE PROCEDURES 
A. Kink motion 

Let us first consider update procedures that are straightforward generalizations of those known in the discrete-time 
worldline algorithm, and work with closed trajectories only. The simplest process involves transformations that do 
not change the number of kinks, but change their types, time positions, and temporal ordering,E3 
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(a\ Q ai (T ai )<5 a2 (r Q2 ) . ..Qa n (T a „) It) — ► (« I Qb 1 ( T b 1 )Qb 2 (n 2 ) ■ ■■Qb n (n n ) It) 



(3.1) 



The number of operators involved in the transformation, their types, and time positions are not constrained, except 
that the two configurations have nonzero weight. Obviously, one could suggest many different realizations of Eq. ( ft .if ), 
and some might work more efficiently than others, depending on the system. Here we describe the procedure called 
"kink motion" ; other procedures have too much in common to be described separately, and allow trivial modifications. 

To move a kink, we first select it at random from the list of existing kinks and decide on the time interval to be 
considered. Suppose that we have chosen a transition described by (Q , t ). We then find kinks of the same type that 
are nearest in time (both to the left and to the right of r Q ), i.e., Q Q or Qj, and consider their times tj < r Q and t 2 > r Q 
as the boundaries of the time "window" transformed by this procedure (in certain configurations at high temperature, 
it may happen that (n, t 2 ) = (0,(3) ). It is allowed to have any number of kinks of different types Q a ^ Q , QJ within 
(ti, t 2 ). Thus the typical initial configuration has the form 



Q (t )... Q an (r a „) 



(3.2) 



(as explained above, one has to consider only those kinks which do not com mut e with Q Q ). 

The second step is to analyze all possible configurations obtained from (3.2) by removing Q from point r and 
inserting it at arbitrary t' £ (ti,t 2 ). We keep the time positions and the chronological ordering of all the other 
operators Q ai , Q ai , ■ ■ ■ , Qa„ untouched. The new position of the selected kink Q D in time is decided according to the 
statistical weight of the final configuration as defined by Eq. (2.3). This is done in complete analogy with the classical 
MC procedure of taking multidimensional integrals. 

The acceptance rate of the kink motion procedure is unity, since the differential measure of the initial configuration 
is zero. In this way, all noncommuting kinks in the Hamiltonian (except kink - antikink pairs, which are dealt with in 
the next subsection) can change places. In dimensions d > 1, the kink motion procedure must be supplemented with a 
"local loop" procedure, which generates small loops in real space, e.g., by replacing Q^i+ gi {T\)Qi+ gi ^i+ gi + g2 (r 2 ) — > 
Qi^i+g2( T 3)Qi+g 2 ^t+ gi + g2 { T i), where gi,g 2 are the nearest neighbor indices. 



B. Creation and annihilation of kink - antikink pairs 

In this subsection we make use of the general theory of Sec. |l] and explain how the elementary procedure of 
creation and annihilation of kink - antikink pairs is organized in practice. An important new principle realized in our 
algorithm is the possibility of selecting different update procedures with certain probabilities (see also Appendix A). 
These probabilities, p a and p c , are at our disposal, and if necessary, can be used to "fine tune" the efficiency of the 
MC process as a whole. The most natural starting point for the update is to address at random some configuration 
fragment. It can be characterized by the kink Q {t ), or by the system state |a(z )) between the two adjacent kinks 
that change this state (in computer memory, all \a(i )} between kinks are assigned labels; the configuration itself is 
described as a linked graph by specifying nearest-neighbor associations (in space and time) between the labels). We 
choose the latter variant and address site labels. Thus the probability of applying an update procedure to a given 
fragment is oc 1/iVjab; where iVi a b is the total number of labels characterizing the initial configuration. By inserting 
(deleting) n extra kinks, we increase (decrease) Ai a b by X)j=i m Qj wnere m Qj gives the number of states changed 



by the kink Qj. Thus the ratio p a /Pc in Eq. (2.7) is proportional to A^ab/C^ab + 2j=i m Qj) when addressing the 
creation of n kinks, and (M a b — Ylj—i TO Qj)/Ai a b when addressing the annihilation procedure. 

To fix the values of p a and p c , we count the number of possible kink - antikink processes that can be applied to 
a given fragment. This number is denoted by N ploc . The simplest choice is then to assign equal weight, l/A proc , 
to all possibilities. For example, if we consider a model with the nearest-neighbor hopping in ID, then there are 
three possibilities for the site state \a(i)): to insert Qi^i+iQi+i^i or Qi^i-iQi-i^i and to delete a pair of kinks 
that change this state to the left and to the right in time, provided they form a kink - antikink pair (i.e., are of the 
Qi±i^iQi—>i±i type). In this case ttlq. — 2 as well, and we finally have 

Pa Ai ab p a A lab - 4 

— = (creation) ; — = (annihilation) . (j-'j) 

Pc Mab +4 p c Alab 

Obviously, in the thermodynamic limit and at low temperature, these ratios are very close to unity. Again, this is 
only a particular example; other choices may prove to be more efficient under certain conditions. 

Once the configuration fragment and update procedure are selected, we proceed along the lines described in Sec. ||. 
Here we would like to comment on the choice of probability density W(f). It would be perfect from the acceptance 
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rate point of view to take W(t) oc R(f). However, this can turn out to be a very expensive procedure. To illustrate 
the point, consider a configuration fragment of length 77. r = t t — 77. Due to the large interaction radius between 
particles, an effective field acting on updated states can change many times during 77 , r . If the number of time slices 
thus induced on the interval (rj,r r ) is N Tlr 1, then complete parametrization of the R(t) function will require 
calculation of the N T[ r (N Tl r + l)/2 partial probabilities, according to the number of ways one can distribute two 
kinks among N n r time subintervals. 

The solution of the problem lies in choosing W(f) = W(E, t), where W(E, f) is an analytic function with the same 
properties as W(t), controlled by a parameter E that is used to minimize the variance of |W(i?, f) — R(t)\. The most 
obvious physical choice of E is the mean field potential acting on the updated states during 77 r from the rest of the 
system 

e-^.r =R(n,T r ) , (3.4) 



p — E(T2— Tl) r T r fT 2 _ 

W{E,f) = - = , /= / dT 2 dne- E ^- T ^ . (3.5) 

One immediately recognizes in W(E, f) the statistics of the kink - antikink pair in the biased two-level systemS, 
which, through the mean-field definition of the bias energy E, most closely approximates the local statistics of kink - 
antikink pairs in a real system. 

The procedures described in the last two subsections represent a direct generalization of local procedures already 
known in the discrete-time worldline method. Their continuous-time versions are, however, only specific realizations 
of a much wider class of possible procedures, thus making the overall CTWL scheme more flexible. 



C. Creation - annihilation, jump, and reconnection procedures for worldline discontinuities 

Up to now, we have considered procedures for working with closed worldlines. These are sufficient to simulate 
quantum statistics in the canonical ensemble and in the M = sector. To overcome this essential drawback, and to 
calculate the Green's function, one usually introduces an extra worldline segment and simulates quantum statistics in 
the presence of two worldline discontinuities at points (ii,T\) and {12, T2)- This process is highly inefficient, because 
one has to probe all degrees of freedom in the configuration (numbering roughly ~ L d j3) to collect statistics for only 
two extra degrees of freedom. In practice, this method was never used to calculate Green's function in large systems, 
e.g., with L d (3 ~ I0 4 . The solution we find for this problem is in considering the two worldline discontinuities to be 
real dynamic variables in the Hamiltonian, which are allowed to move through the configuration both in space and 
time. It turns out that this motion can be arranged to be ergodic, and probes all possible system states. One can 
even completely ignore all the other update procedures, such as moving other kinks and working with kink - antikink 
pairs, probably at the expense of being less efficient, but still remaining accurate, complete, and universal. Below 
we describe the details of update procedures with worldline discontinuities ("worms"), which were first introduced in 
RefM 

We start with the general expression for the Matsubara Green's function (see, e.g., RefrJ) in the interaction picture 



G(i,j,Tx,T 2 ) 



s^Tr 



e-e H °T T ( ai (n)at(T 2 )a) 



(3.6) 



where T T is th e r-ordering operator, which was explicitly written before in defining the Matsubara evolution operator 
fj in Eq. (2.1); £1 is the grand canonical potential. To be specific, we assume here that H Q is diagonal in the site 
representation; in the general case, one might imagine that the index i refers to some parametrization of eigenstates 
of H a . Since we now work in the grand canonical ensemble, the Hamiltonian contains an extra term 



(3.7) 



where [i is the chemical potential. Formally, the only difference between the statistics given by Eq. (2^3) and the 
Green's function (3.6) is that we have two extra kinks, di(ri) and a] (72). Hence one has the possibility of calculating 
the Green's function in a unified process, together with standard thermodynamic averages ("energy," for the sake 
of brevity). To this end, it is necessary just to work in an extended configuration space, where two classes of 
configurations are present: (i) with continuous worldlines, and (ii) with two worldline discontinuities, corresponding 
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while those of the class 



to the kinks aj(n) and aj(r 2 ). (Clearly, configurations of class (i) contribute to "energy,' 
(ii) contribute to the Green's function.) The transitions between the two classes are performed by the processes of 
creation and annihilation of the kinks <Xj(Ti) and at(ra), in accordance with the general balance principles Eq. ( [2.7] ). 
For computational purposes, it is reasonable to redefine the Green's function by a trivial scaling transformation 



77 ai 



na'-, where the constant 77 is adjusted to produce the optimal acceptance (rejection) probability. 



Alternatively, one can arrive at the above scheme by the standard trick of introducing a source to the configuration 
action S (the notation ?/ for the source is chosen deliberately): 



f drV{r) — > f drV{r) + V f ' dr( V *(r) a,(r) + Vi (r) a+(r)) , 
Jo Jo „■ Jo 



(3.8) 



and defining the Green's function as a functional derivative of the generating functional (the partition function with 
the source) 



1 



G(i,j,n,T 2 ) = -— 



S 2 Z 



Z 5th(ti)5t$(t2) 



(3.9) 



The numerical procedure equivalent to the variational derivative in the limit 77 — s- means that only configurations 
with (i) zero and (ii) two worldline discontinuities are included in the statistics. Confining ourselves to just these 
configurations, we do not have to deal any longer with infinitesimally small rj , and can choose rj to be a certain 
finite constant. (This is crucial for any realistic computational process, since 77 — > clearly means that the time of 
accumulation of statistics goes to infinity.) Indeed, a particular value of 77 just defines the relative weights of classes 
(i) and (ii), thus changing the relative norm of the Green's function with respect to "energy" by the known factor 



of 



(Incidentally, one may pay no attention at all to the normalizing statistics for the Green's function, as the 



norm can ultimately be fixed by the condition Q{i, i, t, t + 0) = —density.) 

A typical configuration with two "worms" is shown in Fig. 1 ("live" picture taken from the computer). To update 
it we apply the following transformations: 

• Creation and annihilation of two worldline discontinuities 
We delete a pair, ai(r{)a\{T2) or a\(ri)ai(T2) , when discontinuities happen to meet at the same site i and there 
are no other kinks between them that can change the state of i. The only difference between this and the kink - 
antikink procedure is that now we transform only a single-site state, thus vtlq — I. The annihilation procedure 
addresses the pair of worms, but the creation procedure (which makes sense only when there are no worms) 
addresses the randomly selected configuration fragment label. The ratio of probabilities p a /Pc to address update 
procedures that transform the same configuration fragment back and forth is now 



Pa 
Pc 



= Mab (creation) ; 



Pa 
Pc 



= Mab — 2 (annihilation) 



(3.10) 



This ratio is macroscopically large, which is obviously unpleasant for the computational process. However, we 



have i?(r) 



n 



with the freedom of choosing 77. By setting | 77 | 2 ~ l/(M ab ), where (Mab) is the average 



number of labels in the configuration, we obtain an update procedure that is not based on small parameters (in 
practice, any rou gh est imate like {L d (3) for (Ni^) is sufficient). The rest is done in exactly the same manner as 
described in Sec. IIIB. 



• Jump 

This update procedure is illustrated in Fig. 2. We select one of the worldline discontinuities and suggest shifting 
it in space by inserting an ordinary kink (hopping operator) to the left (in time) of the annihilation operator 
and to the right of the creation operator. As a result, the worm "jumps" to another site. The number of kinks 
changes by one in this procedure, but p a /Pc is unity, because we address it upon the availability of worms, and 
not according to the number of labels. Also, since we are dealing with only one extra kink here, the structure 
of the R(t) function (see Sec. |J) is much simpler, and we choose W(t) = R(t)/ J dri?(r). The integral is over 
the time interval of the updated fragment. The opposite procedure is called an "anti - jump." 

• Reconnection 

Formally, this update procedure, which is shown in Fig. 3, is technically identical to the "jump," but now an 
extra kink is inserted to the right of the annihilation operator and to the left of the creation operator. We still 
distinguish between them, because in the jump procedure the corresponding particle trajectories do not exchange 
places, while they do so in the reconnection update. Figure 3 makes it clear that we have effectively reconnected 
worldline segments of different trajectories. Note that in fermionic systems, any reconnection/antireconnection 
procedure results in a change of the configuration sign. 
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• Shift in time 

The motion of worldline discontinuities in time is essentially the kink motion process (see Sec. [II A ). Su ppose 
that we have decided to shift an annihilation kink Q Q = a^. The only difference from the scheme ( |3.2| ) is in 
the definition of the updated time interval. Its boundaries (ti,T2) now correspond to the time positions of the 
nearest left and right neighbors (kinks) of any type that operate on the same state i. Of course other possibilities 
are allowed as well, if one has some physical arguments in favor of, say, extending the time window farther to 
the next-nearest kink, or a kink of a special type. 

The update procedures thus defined comprise an ergodic stochastic process that operates on the entire configuration 
space of the system. All configurations, including those with nonzero winding numbers and different number of 
particles, are accounted for. Extra particles are inserted/removed from the system when af^Ta) makes a complete 
loop in time (relative to aj(ri)), i.e., when ti — t\ changes by multiples of (3. Winding numbers are introduced when 
j — i changes by multiples of L. The key point of our approach is that each local update makes a contribution to the 
Q(i,r) histogram, except rare cases in which there are no worms in the configuration; these configurations contribute 
to the "diagonal" (or conventional) statistics of closed worldlines. Contrary to the standard calculation, we do not 
adjust all degrees of freedom to the current positions of worms, but rather probe and update the whole configuration 
through their motion. This almost trivial modification results in a factor of (L d j3) acceleration of the scheme! 

It is instructive to draw an analogy between the motion of worldline discontinuities and the loop cluster update 
rules. As is_.easily seen, the basic elements of the single-loop LCU method known as "optional decay" and "forced 
transition"!!!] correspond to a particular evolution of the worldline discontinuities ("optional decay" corresponds to 
the "jump" procedure, and "forced continuation" to the "antireconnection" procedure). A closed loop is obtained 
after annihilating the pair a and a\. Notice, however, that in our scheme (i) not only closed loops, but also all 
intermediate configurations are physically meaningful and are included into the statistics; (ii) nothing is based on the 
special structure of the system Hamiltonian; (iii) the update is always local (it is known that acceptance rates for 
large loops become very small when an external magnetic field in the Z-direction is applied to the Heisenberg system 
(magnetic field is equivalent to a finite chemical potential in bosonic language); this problem is simply absent in local 
schemes). 

The statistics of discontinuities in space and time is given by G{i,r), i.e., it is defined by the Hamiltonian. In 
general, the optimal update scheme depends on the quantity being calculated, and thus one might wish to control 
the statistics of worldline discontinuities "at will." This can be easily achieved by introducing a fictitious spacetime 
dependent potential acting between the "worm" ends, so that their relative positions are now distributed according 
to the function 

Q{i,r) Q(i,r) , 

where Q(i, t) is arbitrary. In this way, one can change the typical size and shape of the loops generated by the "worm" 
algorithm. 

The scope of the present paper is such that we are unable to discuss here many important details concerning the 



practical implementation of our algorithm (optimal triple-linked storage, particular forms of Eqs. (2.8)-(2.9) for each 
subprocess, optimal management of subprocesses, etc.). Readers interested in these issues are encouraged to take 
advantage of our FORTRAN code with comments. The code is written for the ID boson Hubbard modelcJ. 



IV. ILLUSTRATIVE RESULTS 



To demonstrate the advantages of the CTWL algorithm, we have calculated properties of the ID boson Hubbard 
model (Eq. (1.1) with Uj = Uq 5ij) for various coupling constants Uq and particle densities p. 

Comparison with the exact diagonalization results for small systems has demonstrated the lack of any detectable 
systematic error. In particular, for a system with eight lattice sites and six bosons, and on-site repulsion U = 0.5, the 
exact diagonalization result for the ground-state energy is Eq — —10.49209, while long-run Monte Carlo simulations 
yield Eq — —10.4922(2), i.e., a result with relative accuracy better than 10~ 4 . 

It is well known that a commensurate system with p = 1 undergoes a superfluid - Mott-insulator transition of 
the Bcrezinskii - Kosterlitz - Thouless typeao when-the on-site interaction is strong enough (for the most accurate 
estimate of the transition point Uq — 1.645i, see Ref.EZI). In the superfluid phase, including the critical point, one can 
utilize knowledge of the long-wavelength behavior of the system. As explained by Haldaneo, the energy associated 
with extra particles and nonzero winding numbers is quadratic in M and N — N (for simplicity, in what follows 
we count particle numbers from the commensurate value: N — > N — L and N — » N — L). This means that the 
corresponding probability distribution in M and N is a Gaussian, i.e., 
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W[N, M) oc exp 
ex exp 



2/3A s (0) 



M 



P 



(N - N) 2 



2Ln(0) 



(4.1) 



TT\/A S K 



The zero argument of the superfluid stiffness A s and compressibility k denotes values at T = 0. Here K~ 
is the index that controls the asymptotic behavior of the correlation functions, and c is the speed of sound. 

At the critical point, K(L), k(L) and A S (L) are system-size dependent quantities, with K(L — > oo) — * 1/2. Since 
the speed of sound is unrenormalizable in a homogeneous system, it is sufficient to study scaling equations for the 
critical index only. In fact, the solution of the renormalization group (RG) equations for K{L) can be "visualized" 
by considering the logarithmic derivative of the Green's function, since its index is just K/2: 



K(l) = -2 d -^l; l = lnr. 
amr 



(4.2) 



Here we have introduced the variable r 2 
both in space a nd t ime. 



(ct) 2 , which by conformal invariance describes asymptotic decay of Q 



Expressions (4.1) and (4.2) allow for a comprehensive test of the new algorithm. It is also tempting to consider a 
large system right at the quantum critical point and to evaluate its properties under the most unfavorable conditions 
for the standard worldline method. To calculate the critical index and the speed of sound, we considered a ring with 
100 lattice sites and [3 = 100/i. The critical parameters of the Hamiltonian are Uq = 1.645i and /i = 1.94£3. We 
had no problems in accumulating sufficient statistics of winding numbers and N for this system (the corresponding 
calculation is virtually impossible using the standard worldline algorithm). Simple manipulations with the exponents 
in (4.1) result in the following expressions: 



N N ^ rN 



2(N 2 + r N ) 
(3 N 2 



= T ' P N 

L p N 



A s (0) 



L M 2 

9M 



9M 



111 

-In 



\n[W(Q,N)/W(0,0)} 
ln[W{0,-N)/W{0,0)] 
W(0,N)W(0,-N)' 
W 2 (0,0) 

W(0,M)W(0,-M) 



IF 2 (0,0) 



If one is interested in evaluating directly K(0) then 



K(0) 



n\NAl\ 



(4.3) 



(4.4) 



The choice of N and M here is arbitrary, but for numeric reasons, the optimal N and M correspond to values 
where {pn,9m) ~ 1- The advantage of working with nonzero winding numbers in the grand canonical ensemble is 
obvious: in a single MC calculation, one collects all the necessary information about the parameters in the effective 
long- wavelength action, which is very convenient in determining quantum critical points from K — K c . For the 
aforementioned system we found c/t = 2.4(1), and K(l = ln (100)) = 0.47(1). 



One note is in order here. The Gaussian distribution (4.1) implies that the system is in the superfluid phase. In the 
general case one has to define the compressibility as n — dp/dfj,, where by definition p = N / L. The superfluid stiffness 
A s is defined as the coefficient relating persistent current and gauge phase when ip — * 0; this yieldsE3 A s = M 2 L/f3. 

Finally, we used our method to evaluate the Green's function Q(i, r) and to extract the critical index of the 
Berezinskii - Kosterlitz - Thouless transition from its asymptotic behavior; one can then check the consistency of 
all calculations. Since the CTWL simulation yields a two-dimensional histogram for Q(i,r), much more accurate 
results for K(l) are obtained by computing logarithmic derivatives along different directions in the [x, r) plane with 
subsequent angular averaging. The speed of sound, which is necessary for such a calculation, is extracted from the 
asymmetry between x and r in the asymptotic decay of Q. For the Green's function calculation, we considered a 
ring with 450 lattice sites and (3 = 200/i. In Fig. 4, we show the short-range behavior of the Green's function in the 
r-direction, with the characteristic jump at r = 0. The dashed curve is a linear interpolation between the calculated 
points. In Fig. 5, we present the full-scale behavior of Q(x, r) by plotting it as a function of r = (x 2 + (ci) 2 ) 1 / 2 along the 
time and x = ct directions. In accordance with conformal invariance, for large r the two curves are indistinguishable 
to within the statistical errors. The speed of sound obtained from the Green's function is c/t — 2.4(1), and analysis 



of the logarithmic derivative (4.2) yields K(l = ln(100)) = 0.46(2) 
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It is worth noting that our calculations for W(N, M) and Q were performed on a Pentium-90 PC. None of these 
results (e.g., for L > 100) can be obtained by other methods, even with the use of supercomputers. 

The strong on-site disorder at low temperatures is a severe trial for most Monte Carlo schemes. Cluster methods 
suffer from inefficient global Metropolis updates here, while standard canonical-ensemble algorithms suffer from slowing 
down due to one-particle local minima in the effective action (the lowest single-particle states are well localized, and 
probing different configurations requires deep sub-barrier motion) . The unique feature of our "worm" update method 
- the possibility of locally seeding an extra world line at any point in the spacetime continuum - obviates these 
problems. 

To demonstrate the efficiency of our method, we present the results of just one nontrivial calculation - the depen- 
dence of the average particle number on the chemical potential in the Bose glass (BG) phase of the ID disordered 
Hubbard model, Fig. 6. We consider a system with L = 60 sites at j3 = 150, U = it. Disorder is introduced by 
randomly distributing the on-site potential 'Yl,i t i n i between —A and A, with A = 6t. The curve (iV)(/i) allows 
precise determination of the low-energy quasiparticle spectrum of the system to order O.Olt, which is equivalent to the 
calculation of the total system energy to a relative accuracy of order 10~ 4 . The entire plot of (JV)(yn) was obtained in 
a few days of CPU time on a PentiumPRO-200 processor. 

To obtain further evidence of the effectiveness of our method in more familiar problems, the reader is referred to 
a calculation of the superfluid - Bose glass - Mott insulator phase diagram for the 1-D disordered boson Hubbard 
modelEi 



V. CONCLUDING REMARKS 



Although the CTWL algorithm developed here is quite general, some aspects deserve special discussion. What if 
the interaction radius r Q is large? All of the procedures update configuration fragments with the typical duration 
T r — ti ~ 1/t. Since the method is exact, we trace all the kinks within the interaction radius, because they contribute 
to the function R(f). This means that the interval (77, r r ) is further split into N n r ~ r^z ^> 1 subintcrvals (z is the 
coordination number), and each subinterval requires special consideration. If r a is as large as the system size, then 
the whole scheme is in trouble, becoming a "victim of exactness." nnn 

The idea of solving such a problem is demonstrated by the stochastic series expansion methodld'EHij. One might well 
wonder why continuous-time schemes, which contain as an essential ingredient an evaluation of time integrals, work 
as efficiently as SSE, which has all these integrals being evaluated exactly right at the start? Also, why is keeping the 
potential energy U in the r-exponent not at all an advantage if U ~ t? The point is that the MC process is exact 
only for asymptotically long computation times, and there is no reason to calculate anything more precisely than the 
unavoidable statistical error, especially if the corresponding calculation becomes the bottleneck for the whole scheme. 
Evaluating time integrals in CTWL or reproducing exponents by expanding in power series in U ~ t are just two 
cases that illustrate this point. 

Suppose that all particles in the system interact with one another, so that formally r Q = L, but J dr p(r)U(r) = 

F(L) 7^ 00. We divide the interaction Hamiltonian into two parts, HQ(r < r t ) + H-Q(r > r t ), by introducing the 
truncation radius 

/ dr pU(r) = t . (5.1) 

J\r\>r t 

We then write H Q = and combine with V (see Sec. |n|), i.e., the lon g-ra nge part of the interaction Hamiltonian 
is now considered to consist of diagonal kinks. Because of the definition ( |5.1| ), the total number of kinks within the 
time interval ~ 1/t remains finite and independent of system size. 

The case of a divergent integral J dr p(f)U{r) — F(L) — > 00 is more subtle, since the number of diagonal kinks 
within the time interval ~ 1/t, given by F(L), is now large (if F(L) is a logarithmic function of L, we do not regard 
this problem as serious). On the other hand, for long-range interactions the so called "mean- field approximation" 
becomes more accurate. Since the mean- field potential is easy to account for analytically (and numerically), one now 
has to deal with fluctuations, and these quite often satisfy the condition 



dr (p(f) — p) U (r) 



SF(L) ^ 00 . (5.2) 



In Appendix A we explain how to organize the Monte Carlo process using the mean-field approximation for the 
configuration weight. The net result is that even for long-range potentials, the calculation time can remain independent 
of system size. 
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In this paper we have concentrated on the Green's function calculation by restricting the number of worldlinc 
discontinuities to or 2. Of course the scheme can be trivially extended to include the case with a larger number of 
discontinuities, if one is interested in the two- or n-particle Green's function or n-point vertex. More generally, our 
scheme makes it possible to work with Hamiltonians that do not conserve the number of particles, i.e., when there 
are sources with finite strength in the bare Hamiltonian. 

Although in this paper we consider a system with discrete Hilbert space in detail, the principles of update in 
continuous time developed here are much more general. Mathematically, we construct an exact method (in the 
statistical limit) of averaging over a distribution represented as a series of integrals with an ever-increasing number 
of variables, but with essential similarity among the terms of the series, allowing their local comparison (weighting). 
We may call such structures integrals with a variable number of variables - VNV integrals. Physically, we sum a 
perturbative expansion in the interaction picture for some observable of a large but essentially finite-size system. 
(For a system with discrete Hilbert space, the only continuous variables in this expansion are the times of virtual 
transitions.) But perturbative expansions for continuous systems also have the structure of VNV integrals, with 
additional integrations over some continuous variables. Thus (apart from the fact that for spatially continuous 
systems one cannot expand the kinetic part of the Hamiltonian and must use the potential energy as a perturbation), 
there is no qualitative difference between perturbative exp ansio ns fo r c ontin uous and discrete systems. The general 
method of evaluating VNV integrals is given by Eqs. ( [2.5| ) and (2/7) - ( [2.9|) , where the vector f now stands for any 
set of continuous variables, and the function R(t) is defined straightforwardly, given the particular form of the series. 
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APPENDIX A: LONG-RANGE POTENTIALS 



Suppose that we are dealing with the case F(L) — > oo, but finite 5F(L). The idea is to organize the Monte 
Carlo process in such a way, that in most updates we simply ignore fluctuations, and account for distant particles 
by replacing them with a homogeneous density distribution. Obviously, for the scheme to remain accurate, in some 
updates we have to consider deviations from the mean-field distribution. The goal is to address the procedure dealing 
with distant fluctuations with the small probability which is at least inversely proportional to the number of operations 
in this procedure. 

Consider again the balance equation for the given pair of subprocesses, but now including the possibility of com- 
pleting the same update procedure in a number of ways: 

A oPc W(r) J2 7 0) P^ c (r)df - dA n (r) Pa £ 7 (j) PgLif) = . (Al) 

j=0 j=0 

Here is the probability of using the j-th version of the update procedure. We require X)j=o 7 = ^ ana - 
7o 71 . . . ^> 7j, . We also assume that the proced ure j» corresponds to the exact treatment of all fluctuations. 



Other quantities have exactly the same meaning as in (2.5). The self-balance condition now reads (compare Eq. (|2.7|)) 



W(f) J2 7 (j) P£l(7) = R(r) £ 7 0) P^(f) . (A2) 
j=o j=0 

To satisfy ( |A2| ) we suggest the following scheme. Let R^\f) be the distribution corresponding to the exact 
treatment of fluctuations up to the distance r^' with <C <C • ■ • <C r^*^ = L, and the mean-field treatment of 
more distant (r > r^') particles. We can write then 

R U) = R W + SR W + 6R U) ( RV-)=R. (A3) 
If SF is finite and is sufficiently large, then all SR^ are small. We then choose 7' ' w 1 and 
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P^lir) = ( R ( y)/W{r) , * < WW , (A4) 

acc ; 1 , otherwise , v ' 



p(o) = f ^(r)/i?(°)(f) , if i?(°)(f) > W(f) , 

romV 7 I 1 , otherwise , v ' 



and solve the self-balance condition deductively by requiring 

k 



W{f) J2 1 U) P^l(r) = R {k) (r) £ 7 0) ^(r) , (A6) 



or equivalcntly 

fe-i 



7 W [ W P W _ pW P W] = 5 pW ^ 7 W) pW . ( A7) 



The final answer can be written 



[JiZW^Io 1 7 (J) Pr ( e^]/[7 (fe) W(f)l , if 5RW > , 
PiS = { " (A8) 

, otherwise , 



f -[^Eto 1 P£2j/[7 (fc) a (fc) (7)] , if ^ (fe) < , 
= (A9) 
[ , otherwise 

Since all SR^ are assumed to be small, it is possible to keep 7^' <C 1 (for k = 1, 2, . . . j*), but large enough to avoid 
situations with Pici > 1 or > 1. 

FIGURE 1. 

A typical 8-site configuration with two worldline discontinuities marked by filled circles. The width of the solid line 
is proportional to the site occupation number, and dashed lines are empty sites. 

FIGURE 2. 

Jump procedure for the annihilation operator, a) Initial configuration fragment; b) suggested variation (in the 
antijump procedure, (b) is the initial configuration and (a) is the suggested variation). 

FIGURE 3. 

Reconnection procedure for the annihilation operator, a) Initial configuration fragment; b) suggested variation (in 
the antireconnection procedure, (b) is the initial configuration and (a) is the suggested variation). 

FIGURE 4. 

Short-time behavior of the Green's function 0(0, r) of the commensurate 1-D Hubbard model at the quantum critical 
point. 

FIGURE 5. 

Long-range behavior of Q(i,r), demonstrating conformal invariance. 
FIGURE 6. 

Number of particles vs. chemical potential in a large Bose-glass cluster at macroscopically low temperature. 
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